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Introduction 


To  facilitate  the  OTHR  coordinate  registration  process,  bistatic  backscatter 
sounders  are  routinely  operated  along  with  QVI  sounders  at  the  two  operational  facilities 
in  Virginia  and  Texas  [Headrick,  1990;  Headrick  and  Skolnik,  1974].  These  sounders 
have  azimuthal  receive  antenna  beam  steering  capability,  so  that  a  BI  can  be  measured  for 
a  set  of  beam  directions  within  the  azimuthal  coverage  sector  of  each  backscatter  sounder. 
Here,  we  combine  each  of  the  eight  individual  WSBIs  obtained  at  different  azimuths  to 
form  an  integrated  view  of  the  ionospheric  density  over  the  OTHR  coverage  sector  (64°- 
80°). 

Each  WSBI  typically  has  a  distinct  leading  edge  which  is,  to  first  order, 
independent  of  the  scattering  properties  of  Earth’s  surface  and,  in  the  case  of  a  narrow 
antenna  beam,  is  determined  entirely  by  the  three-dimensional  distribution  of  electron 
density  existing  in  the  ionosphere.  The  theoretical  leading  edge  for  a  specific  ionospheric 
model  may  be  found  to  high  accuracy  by  using  ray  tracing,  as  discussed  in  Section  2. 

The  leading  edge  of  the  WSBI  contains  information  about  ionospheric  regions 
located  thousands  of  kilometers  away  from  the  sounder  and  the  inverse  problem  of 
backscatter  sounding,  as  we  address  it  here,  is  the  problem  of  extraction  of  the  plasma 
density  distribution  in  the  ionosphere  from  the  measured  leading  edge.  This  inverse 
problem  is  rather  cumbersome,  because  even  the  solution  of  the  direct  problem  (the 
relationship  between  the  plasma  density  distribution  and  the  leading  edge)  can  not  be 
expressed  explicitly.  A  rigorous  approach  that  addresses  the  numerical  solution  of  this 
problem  was  developed  recently  by  Fridman  and  Fridman  [1994].  This  approach  is 
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based  on  the  Newton-Kontorovich  method  for  the  solution  of  nonlinear  operator 
equations  and  on  Tikhonov’s  regularization  method  for  ill-posed  problems.  Fridman  and 
Fridman  [1994]  considered  the  inverse  problem  for  a  single  WSBI  and  reconstructed  the 
two-dimensional  electron  density  distribution  in  a  vertical  plane  aligned  in  the  direction 
of  sounding.  The  technique  was  found  to  accurately  reconstruct  the  ionosphere  in 
numerous  tests,  both  against  simulated  data  as  well  as  multi-station  experimental  data 
{Fridman  et  al,  1994]. 

In  this  project,  the  approach  of  Fridman  and  Fridman  [1994]  is  extended  to 
develop  a  method  for  the  reconstruction  of  the  three-dimensional  ionosphere  from  the 
OTHR  WSBI  data.  A  regularized  formulation  of  the  WSBI  inversion  problem  and 
specific  numerical  methods  for  this  nonlinear  inverse  problem  are  described  in  Section  1 . 
The  key  element  for  the  solution  of  nonlinear  problem  is  Tikhonov’s  method  for  two- 
dimensional  linear  problems  developed  in  Section  2.  Section  3  presents  several  examples 
of  inversions  for  data  collected  by  the  OTHR  located  near  Chesapeake,  VA.  In 
conclusion,  the  WSBI  inversion  technique  is  discussed  as  a  new  method  for  ionospheric 
diagnostics. 


2 


1.  Theory  of  WSBI  Inversion 


1.1  Starting  Principles 

The  measured  WSBI  leading  edge  is  a  function  of  two  variables  g(/,(p)  because 
it  represents  the  dependence  of  the  minimal  group  delay  on  frequency  /  and  azimuth  (p . 
The  theoretical  value  of  this  function  g(/,(p),  for  any  given  distribution  of  electron 
density  n{xi,X2,X2,)  specified  in  a  system  of  coordinates  x^,X2,Xt,,  may  be  calculated 
by  munerical  solution  of  the  ray  tracing  equations  as  follows:  for  each  frequency  and  each 
azimuthal  beam,  a  fan  of  oblique  rays  exiting  in  the  central  azimuthal  direction  of  the 
beam  at  different  elevation  angles  is  considered,  and  the  ray  that  exhibits  the  minimal 
group  delay  for  propagation  from  the  sounder  to  a  point  of  first  contact  with  the  Earth  and 
back  to  the  sounder  is  derived.  This  minimal  group  delay  and  the  current  value  of 
operating  frequency  represent  two  coordinates  of  a  point  in  the  leading  edge  plot.  Every 
point  of  the  leading  edge  is  obtained  by  applying  the  above  procedure  to  the  full  set  of 
operating  frequencies  (the  frequency  range  is  5-28  MHz  for  the  Virginia  OTHR). 

Formally,  this  process  establishes  a  nonlinear  operator  G„that  transforms  n{xi,X2,x-i) 
into  g(/,(p).  It  can  be  represented  symbolically  as 

g  =  \n.  (1) 

Essentially,  WSBI  inversion  may  be  viewed  as  resolving  this  equation  with  respect  to  n. 
The  equation  is  obviously  imderdetermined  because  the  known  fimction  g  is  a  function 
of  two  variables,  whereas  the  function  to  be  found  depends  on  three  variables.  In  order  to 
avoid  this  inconsistency,  we  parameterize  n  by  a  function  of  two  variables  that  is  denoted 
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as  u(5i,52).  The  physical  meaning  of  u  and  of  the  variables  5'i,j'2  depends  on  the 
selected  form  of  the  parametrization.  The  most  general  parametrization  may  be  formally 
represented  as 

n  =  Nu  (2) 

where  iV  is  a  known  operator  (generally  nonlinear)  that  transforms  u  into  n.  The  method 
described  below  is  designed  to  work  with  an  arbitrary  parametrization,  however  in  all  the 
examples  presented  in  this  report,  we  use  the  following  simple  form  of  parameterization 

n(x, ,  X2 ,  X3 )  =  no  )[1  +  =  ^2 )]  >  (3) 

where  no  (X3 )  is  the  vertical  profile  of  electron  density  as  measured  by  the  QVI  sounder 
at  the  OTHR  site.  Here  and  below  we  require  that  coordinates  Xi,X2  specify  the 
geographic  position  of  a  point  and  the  altitude  is  determined  by  X3 . 

After  the  parametrization  is  introduced  we  have  g  =  G„M,  where  =  G^N 
defines  the  group  path  operator.  The  inverse  problem  is  the  task  of  resolving  the 
functional  equation 

g  =  GuU,  (4) 

where  u  is  the  unknown  function  in  this  equation. 

Inverse  problems  such  as  the  one  considered  here  are  usually  unstable  [Tikhonov 
and  Arsenin,  1977],  meaning  that  even  the  smallest  error  of  measurement  in  g  causes 
considerable  deviation  of  the  solution  u  from  the  true  solution.  In  order  to  stabilize  the 
solution  we  have  incorporated  Tikhonov’s  method  of  regularization  [Tikhonov  and 
Arsenin,  1977]  in  this  work.  The  essence  of  this  method,  as  applied  to  the  WSBI 
inversion  problem  may  be  formulated  as  follows:  All  possible  functions  u  that  satisfy 
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equation  (4)  within  the  accuracy  of  measurement  of  the  leading  edge  g  are  considered, 
and  among  these  functions  the  smoothest  one  is  chosen  as  the  optimal  solution. 

In  order  to  implement  this  approach,  a  quantitative  measure  of  the  smoothness  of 
u  must  be  introduced.  This  measure  is  called  the  stabilizing  functional  Q[m]  and  is 
defined  by 


»  * 

Q[w]  = 


ff 

2  1 

2' 

IJ 

U  +9il 

iwdsids2  , 


(5) 


where  qi(si,S2),  q2(si,S2),  and  w(5i,52)  are  arbitrary  positive  weight  functions  that 
are  chosen  in  accordance  \vith  the  symmetry  of  the  system  of  coordinates  and  the  typical 
scales  of  the  task  in  consideration.  Let  6  be  an  a  priori  root-mean-square  error  of 
measurement  of  the  WSBI  leading  edge,  then  the  regularized  solution  of  the  inverse 
problem  obeys  the  following  relationships: 

^max  /in  ax 

J  j(g  ~  d(pd/  <5  (tp  max  ~  9  min  X^max  — /min ) 

/itiin 

Q[w]  ->  min  (^) 

where  cp^iax  =  9  min  and  are  respectively,  the  upper  and  lower  bounds  of  the 

azimuthal  direction  and  frequency  of  the  transmitted  sounder  beam  pattern.  Thus,  the 
functional  Q[m]  is  minimized  on  all  u  that  satisfy  the  inequality  (6). 

A  procedure  for  the  numerical  solution  of  the  problem  stated  by  (6)  and  (7)  has 
been  developed  as  a  generalization  to  three  dimensions  of  the  technique  for  BI  inversion 
in  two  dimensions  described  in  Fridman  and  Fridman  [1994].  It  is  an  iterative  process, 
where  each  iteration  executes  the  following  three  steps: 
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1 .  Ray  tracing  synthesis  of  the  WSBI  leading  edge  for  a  given  ionospheric  model, 

2.  Linearization  of  the  operator  (based  on  the  solution  of  the  linearized  ray  tracing 
equations), 

3.  Calculation  of  corrections  to  the  ionospheric  model  by  solving  the  problem  stated  by 
(4)  and  (5)  with  replaced  by  its  linearized  approximation.  The  solution  is  attained 
by  applying  the  regularization  method  for  linear  problems,  a  method  developed 
originally  by  Tikhonov  and  Arsenin  [1977]. 

The  corrected  model  obtained  from  step  3  is  used  in  the  next  iteration  for  step  1.  The 
iterations  are  continued  until  the  model  converges. 

Input  data  for  the  inversion  proeedure  include:  (i)  a  table  of  WSBI  leading  edge 
i  (ii)  the  apriori  root-mean-square  error  5  of  the  leading  edge  measurement;  (iii) 
the  overhead  vertical  profile  n^ih)  obtained  from  the  QVI  data  {h  is  altitude).  The  output 
is  a  three-dimensional  distribution  of  electron  density  for  the  area  covered  by  the  sounder. 
The  maximum  height  attained  by  this  inversion  is  always  below  the  F2  region 
ionospheric  maximum,  because  all  rays  forming  the  leading  edge  are  refracted  below  the 
peak  of  electron  density. 

In  this  report  we  will  describe  the  most  important  developments  of  the  inversion 
technique  achieved  during  the  course  of  this  project.  In  1.2  we  will  describe  the  theory 
and  algorithm  for  calculation  of  the  linearized  group  path  operator.  Then,  in  1 .3,  our 
realization  of  the  Tikhonov’s  regularization  method  for  the  nonlinear  problem  in 
consideration  is  presented. 
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1.2  Theory  and  Algorithm  for  Calculation  of  the  Linearized  Group  Path  Operator 


Let  us  denote  as  P(P,(p,/)  the  group  path  from  the  transmitter  to  the  point  of  the  first 
contact  with  the  Earth  for  a  ray  exiting  the  transmitter  at  the  elevation  angle  P  (measured 
from  the  vertical  direction),  azimuth  cp ,  and  operating  frequency /.  Available  ray-tracing 
codes  are  able  to  perform  accurate  calculations  of  this  function.  The  leading  edge 
g(/,(p)  is  related  to  the  function  P  as  follows: 

g(/,9)=  min  P(P,(p,/)  (8) 

0<P<7t/2 

The  elevation  angle  at  which  the  group  path  P  is  a  minimum  we  will  denote  as  P^  (9,/)  • 
So  that 

?(/,<(>)  (9) 

Consider  variations  of  functions  P ,  P;„ ,  and  g  associated  with  an  infinitesimal  variation 
hn  of  electron  density  in  the  ionosphere.  That  is  suppose  that 

n(xi ,  X2 ,  X3  )=  riQ  (xi ,  X2 ,  X3  )+  6«(xi ,  X2 ,  X3  ) 

p(p,cp,/)=  Po(p,(p,/)+6P(p,(p,/) 

g(/,9)=  go(/=9)+Sg(/,(p) 

Variation  of  the  leading  edge  5g  may  be  expressed  in  terms  of  6P  and  5p^ .  From  (9) 

we  have:  6g(/,9)=  6P(p^o>9>/)+  5P„((p,/)5Po(P;„o>9,/)/5P  •  The  second  term  in 
the  right-hand  side  is  equal  to  zero  because  of  (8),  so  that 
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5^(/,9)=SP(P;„o.9,/)  (10) 

In  order  to  proceed  further  we  will  need  to  use  a  concrete  system  of  ray  tracing 


equations.  However  the  results  obtained  may  be  easily  generalized  for  any  arbitrary  form 
of  ray-tracing  equations. 

To  simplify  the  procedure  and  to  reduce  the  computation  time,  we  neglect 
azimuthal  displacements  acquired  by  the  rays  during  their  propagation  in  the  ionosphere. 
This  approximation  is  justified  by  the  fact  that  in  the  first  approximation  these 
displacements  do  not  affect  the  leading  edge.  Effects  of  the  geomagnetic  field  will  be 
also  neglected.  With  these  simplifications  the  ray  tracing  equations  may  be  represented 
in  the  form: 


^=F,{x.y^.y^) 

where 


(11) 


(12) 

(13) 

(14) 
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>’]  =  10g(l  + /z  /  i?£;  )  (15) 

y2=4^cos^  (16) 

;;3  =  P'  (17) 

z  =  R^\og(l  +  h/  (18) 

s  =(l  +  /z/i?£)^s  (19) 


8  =  1 


f 


(20) 


/.= 


ne 


Tzm 


(21) 


P'  is  the  current  value  of  the  group  path,  and  x  is  the  range  along  the  earth  surface 
measured  in  the  direction  of  the  exit  azimuth. 

Let  us  write  down  linearized  equations  (1 1)  by  substituting  in  (1 1)-(20)  the 
following  expressions 

=.^10  +6yi 
yi  -  yio  +  6^2 
y3  =  >^30  +  63/3 

/p  =/p0+6/p 

and  implying  that  Jio  ,:>^20 = ^  ^nd  fpQ  satisfy  to  (1 1).  Asa  result  the  following  system 


of  equations  may  be  obtained: 


d 

dx 


5y, 

ail 

"12 

o' 

5>'i 

b\ 

0’ 

63^2 

+ 

^21 

^22 

0 

53^2 

= 

h 

¥p  + 

C2 

6^3  _ 

.«31 

^32 

0_ 

.6^3  _ 

P'i, 

0 

dh 


(22) 


where 
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.  N _ 3^20 _ 

2^/2  0^  So. 

2(^0  -  >’20  J 


a]2(x)  =  — 


l^O  “  >20 , 


^  _ ^ _ g'sQ  1  faeoV 


«22W  =  — 


>20 


■(so  “>20  J 


^2  dz 


«3lW  =  - 


exp(2yiQ/RE)(  4 


1  der 


2^20  ->20  So -3^20 


exp(23^io/^£)v20 
"32W  =  — 7: - 2  V/2  - 

i,So  -  >20  / 

,  >2oexp(2;^io /i?^) 


b,=- 


Qxp(2yiQ/ Re) 


^E  Vso  “>20  2(so  -  >20  J 


1  d  ^ 

,.2  V/2  So 


exp(2>^io/i?£) 

2/^(So-y|or 


ex  p(3>io  ^  ^e) 


For  our  purpose  equation  (22)  should  be  considered  with  zero  initial  conditions.  The 


nearest  goal  is  to  find  the  relationship  between  5g  and  5/^  that  follows  from  (22). 
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Suppose  that  the  value  of  ^  unperturbed  ionosphere  is  found.  Then  we  take  the 

solution  that  corresponds  to  the  exit  angle  the  zeroth  order  solution  to  be  used  in 

(22).  Suppose  that  equation  (22)  is  solved  for  a  given  5/^ ,  then  according  to  (10) 


5g  —  §P(P^0 Ts "t' ^  TsoC^/wo)"*" 
5Ti(^mo)  d 


T3o(^/«o)+ST3(^/«o) 


So  that 


5g  =  6P(p^o)=  ST3(^mo)- 


•^30 

^,«oAy2(->^mo) 

^10 

(33) 


Here  x„j  denotes  the  value  of  x  for  the  first  landing  point  of  the  ray,  and  p  is  the 
elevation  angle  at  arrival  in  the  landing  point  ’  -^10  and  F30  are  correspondingly  F]  and 
Fj  for  the  unperturbed  ionosphere. 

Now  it  is  convenient  to  introduce  Green’s  function  for  equation  (22)  (or  unit 
impulse  response  function).  Namely  we  introduce  a  3  by  3  matrix  R(x,x')  that  satisfies 
the  equation: 


-^R  +  aR  =  I6(x  -  x') 
ax 


with  zero  initial  conditions  (R(jc,x')  =0  at  x  <  x' ).  Here  I  is  a  unit  3  by  3  matrix,  and 


a,i 

^12 

a2\ 

<322 

03, 

^32 
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The  variation  of  the  leading  edge  may  be  written  as 


=  I  ]  ')5/p  [^1  (x X2  (x X2  (x ')]+  G "(x 0^  [xj  (x '),  X2  (x Xj  (x  ')]|^& 


where  the  functions  Xi(x'),X2(x'),X3(x')  are  specifying  the  ray’s  trajectory,  and 


n=\  I  ^10P/«0AT2(^m0jJ 

V,/  A  \  ’  n  /_  A  ^30  ^mO  ATiC^/wo)  r,  /_  A  / 

'  (x  )-  2^  I  ^3«  V^mO  >  X  )-  p_  -  7^  v^wO  >x:  )k„  (x ' 

„=1  [  ^10P;«oAT2(5/n0jJ 


^32(^/nO’^  )  COS^  -^12(5w0»^  )  ^2('^  ) 


The  unit  impulse  response  may  be  expressed  in  terms  of  three  linearly  independent 
solutions  of  (22)  with  6/^=0.  We  obtain  such  solutions  ej(x),e2(x),e3(x)  by  solving 
uniform  equation  (22)  for  three  instances  of  initial  conditions: 


e,(0)=  0  ;e2(0)=  1  ;e3(0)=  0  . 

_oJ  [oj  [l 

Then, 

R(x,x')=  E(x)E''(x')0(x-x') 

Here  matrix  E  consists  of  columns  ej(x),e2(x),e3(x) : 
E(x)=[e,(x),e2(x),e3(x)]; 

and  e(x)  is  the  unit  step  function  (0(x)  =  0  at  x  <  0  and  0(x)  =  1  at  x  >  0). 
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Calculation  of  the  inverse  matrix  is  simplified  by  the  fact  that 


'o’ 

1 - 

0 

1 _ 

63 (x) = 

0 

so  that  E(x)  = 

ei2W  ^22  W  0 

1 

_^u{x)  £23  W  1. 

E-'(x)  = - 

^22^11  “^21^12 


^22 

~^\1 

^22^13  +^21^23 


-£21  0 

£„  0 

^13^21  ~  ^23^11  1. 


(36) 


The  components  of  R  involved  in  the  expression  (34)  for  leading  edge  response  are: 


^12(^m0>^')  = 


gll(^>«o)g22(^')-g2l(^mo)gl2(^0 

e22(^')eil(^')-^2l(^')ei2(^') 


-  gl  1  (^/nO  )g2 1  0  +  ^21  (^mO >^1 1  (^') 

^22  (x  ')ei  i(x')-e2i  (x')ei2  (x') 


(37) 

(38) 


i?j3(x:,x')=0  (39) 

/_  gl3  (Xmo)^22  (x')  -  g23  (^mo)gl2  (^0  "  g22  (^0^13  (^')  +  ^21  (^0^23  (^0 

)  e22(x')eu{x')-e2i(x')e^2ix') 


(40) 


-  ^13 (^mo)g21  jx')  +  g23(^mo)gl  1  (^')  +  ^21  (^0^13  (^')  "  ^1 1  (^'>^23 (^') 
^22  (^0^1  l(jc')  -  ^21  ix')ei2  (X') 


(41) 

(42) 


Let  us  approximate  integrals  and  derivatives  in  (34)  by  finite  sums  and  differences. 
Denoting  as  W/  the  coefficients  the  integration  formula  we  obtain: 
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\G'{x'i)bfp  [xj  (x;)x2  G"(x;)^5/ ^  x,  (x;)x2  ^/z(x/)j|w/  = 


X  W)5/p  [^1  (^/)^2  (^/)K^/)]+ 
/ 


Furthermore,  by  invoking  interpolation  the  above  expression  may  be 
approximately  represented  as; 


=  X  ^)’  (45) 

ij,k 

where  integer  variables  i,  j,  and  k  denote  the  nodes  of  a  3-D  spatial  grid  in  which  the  3-D 
ionosphere  is  specified.  Using  triple  linear  interpolation  we  obtain  the  following 
expressions  for  nonzero  elements  of  the  matrix  G: 

^ijk  =  X  0 "  p\  (oX^  “  Pi  (oX/ 

/:/o(/)=i&yo(/)=7&*o(0=A:‘-  ^-1 

+  ;7i(/)(i-/?2(0X/ 

+  X  ^  ~  P^  — ~7r  C  “  Pi  (^))p2  (Ow/ 

/:/o(/)=/&yo(0+W&*o(0='tL  "*+l 

+  X  ^/^3(o+g/'— ^7 — (^-Piio^-Pii^yi 

Uo{l)=i&Ml)==J&koil)+l=k^  ^  ^-1-1 

+  X  ^^'^3  (0  +  ^iT~\ —  0  “  P\  (0);^2  (J>i 

/:/o(/)=/&7o(0+W&*o(0+1=*L  *  ^-’-1 
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4- 


E 


|■.ii){l)+\=i8^jQ(l)=jSLkQ{l)+\=k 


Glp^{l)  +  G\' 


h  ~  ^k-\  J 


Pi0)(^-P2(i)yi 


+ 


E  [G;(i-P3(0)-Gr^_, 


l:io{l)+l=i&joil)+\=j&W=k 


k+l  ~  '^k 


P\{l)P2Q>l 


E 

/:/o(/)+l=!&7o(0+l=y&^o(0+l=* 


Glp,{l)  +  Gi'- 


1 


K  -  ^k-\ 

The  following  functions  has  been  introduced  above: 


P\{l)P2il>i 


/:o(/):xi(/)  =>io  =^min 

xi(/)>xj'”"’'  =^Zo  =  W  -1 

xj'""’  <  Xi(/)  <  xj'"“  =>  io  =  max/o:x'“  <  Xi(/) 


(46) 


(47) 


Xi(/)-x;°^^^ 

„'o(0+l 

Xi  -Xj 


(48) 


The  pairs  of  functions  joil),P2(0  and  A:o(/),;?2(0  are  defined  in  a  similar  way  for 

spatial  variables  X2  and  X3  = /z ,  correspondingly.  Notations  x|,X2,X3  = /z^  are 
introduced  for  nodes  of  the  variables  Xj ,  X2 ,  X3 . 

When  realizing  these  formula  it  is  convenient  to  have  the  variables  Xj  and  x'  to 
be  the  same.  Then  each  sum  in  (46)  contains  only  one  element  and  the  second  sum  is 
equal  to  zero.  This  also  permits  an  economical  storage  of  the  array  G.  Its  nonzero 
elements  may  be  specified  as  G’zotOJolO+i.^otO  ’  ^fo(OJo(0,^o(0+i’ 

^'o(0  /o(0+i,'to(/)+i  stored  in  four  vectors  of  length  z  max  *min  ^  • 
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Now  consider  the  parametrization  of  by  a  funetion  of  two  variables  u{s,f). 
Linearized  relationship  (2)  may  be  represented  in  the  form: 

6/^  (x^ ,  X2 ,  X3 )  =  j*j*0(xi ,  X2 ,  X3 , 5,  t)?)u{s,  t)dsdt  (49) 

We  approximate  the  above  integral  by  a  finite  sum  with  weight  coefficients  : 

EE  k,  )6u(m,n)r^„ .  (50) 

m  n 

Example.  In  the  case  of  parametrization  (3) 

0(X,  ,X2,X2,S,t)  =  fpQ  (/2)5(x,  -  s)3(x2  - 1) 

and  the  finite  difference  approximation  of  the  above  expression  is 

^(.hjik,S^,tn)  ~  f im^ Jn /^ij 
Using  (50)  in  (45)  obtain: 

6g  =  EE  K'(m,n)r„„du(m,n)  (51) 

m  n 

where 

K'{m,n)  =  (52) 

*  7  k 

Note  that  (51)  is  a  finite  sum  approximation  for  the  integral 

5g  =  j*JjC'(j',r)6w(5,r)<^'S'<^^  (53) 

It  should  be  remembered  that  the  function  K'  depends  also  on  frequency  and  the 
azimuthal  direction  of  sounding. 
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Relationships  (35),  (37)-(42),  (46)-(48),  (52),  (51)  determine  the  desired  linearized 
group  path  operator. 

1.3  WSBI  Inversion  Using  Newton-Kontorovich  Method  and  Tikhonov’s 

Regularization  Method 

Input  data  for  the  inversion  procedure  include:  (i)  a  table  of  WSBI  leading  edge 
g(/,(p) ;  (ii)  the  apriori  root-mean-square  error  6  of  the  leading  edge  measurement;  (iii) 
the  overhead  vertical  profile  Wq  (h)  obtained  from  the  QVI  data  {h  is  altitude),  and  (iv) 
the  apriori  root-mean-square  error  p  of  the  vertical  profile. 

Our  ionospheric  model  is  specified  by  an  unknovm  function  of  two  variables 

u(a'i  ,5'2 )  and  a  known  operator  N : 

«(xi,X2,  hy  Nu  (54) 

The  function  must  be  in  agreement  with  the  measured  WSBIs.  It  also  must  be  in 
agreement  with  «q  (Jt)  ,  that  is  there  should  be 


at  the  QVI  location  specified  by  coordinates  Xio,X2o  •  We  imply  that  the  operator  N  is 
consistent  with  (55).  Thus,  the  functional  equations  to  be  solved  are: 


GuU  =  g 

(56) 

Nsu  =  riQ 

(57) 

Here  we  have  introduced  the  operator  N ^  that  produces  the  vertical  profile  at  the  QVI 
location: 
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(58) 


N^u  =  M/L  ^ 

O  '-*^1 --*10  >^2-^20 

We  are  looking  for  a  regularized  approximate  solution  of  (56),  (57).  We  define  the 
regularized  solution  according  to  the  following  relationships. 

II  ^  l|2  II  ''  l|2 

(1-p) - +  ^ - <1  (59) 


Q[w]  ^  min  (60) 

The  norm  symbol  \\\\  is  used  here  for  rms  value  of  a  function: 


Expression  (59)  indicates  that  we  consider  (56)  and  (57)  as  one  equation.  Parameter  is 
a  positive  constant  smaller  than  1,  it  is  introduced  in  order  to  be  able  to  properly  weight 
these  equations.  Our  typical  choice  of  this  parameter  is  p  =  l/Z',  where  L'  is  the  total 
number  of  points  in  all  WSBIs. 

The  functional  equations  (56)  and  (57)  are  generally  nonlinear.  We  solve  this 
system  iteratively.  In  each  iteration  operators  are  approximated  by  their 

linearizations.  This  approach  to  solution  of  nonlinear  operator  equations  is  known  as 
Newton-Kontorovich  method  [Kolmogorov  and  Fomin,  1975]. 


18 


Linearization  of  the  operator  N  is  specified  by  a  known  matrix  O  in  accordance 
with  (50).  Without  loosing  the  generality  we  may  accept  that  u=0  corresponds  to  the 
uniform  ionosphere  with  the  profile  equal  to  «o(^)  everywhere  and  that 
x,o  =0 

X20  =  0  ■ 

Major  steps  of  the  iterative  solution. 

1 .  The  starting  approximation  to  the  3-d  model  of  the  ionosphere  is  formed  by  assuming 

2.  Ray  tracing  is  used  to  calculate  the  theoretical  leading  edge  g(/,cp)  and  the  optimum 

elevation  angle  (cp,/)  functions  for  all  values  of  frequency  and  azimuth  at  which  the 
experimental  leading  edge  points  are  specified. 

3.  The  rms  discrepancy  between  the  experimental  and  theoretical  leading  edges 


is  calculated.  If  the  condition  (1  -  p)|g  -  g||  ^  6  is  satisfied  then  the  ionosphere  specified 
by  (54)  gives  the  optimum  solution  of  the  inverse  problem.  Otherwise  the  iterative 
process  described  below  is  started  with 
Sg  =  g  -  g 
briQ  (h)  =  0 . 

4.  The  linearized  group  path  operator  K’  (defined  in  (52))  is  calculated  as  described  in 
1 .2  using  P;„(cp,/ )  and  n(xi,X2,h^  found  on  the  previous  step. 
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5.  The  system  of  linear  integral  equations 

(/,  9=  =  6g(/,(p) 

jj®(0  ,0,h,s,ty>u(s,t)dsdt  =5«o(/z) 

is  solved  numerically  using  Tikhonov’s  regularization  method  described  in  Chapter  2. 
When  performing  this  solution  the  system  is  treated  as  a  single  integral  equation  of  the 
form 


II 


K(l,s,t)6u(s,t)dsdt  =  fi(l) 


with  the  discrete  representations  of  the  kernel  and  the  RHS  taken  as: 


K(l,i,J)  = 


K'(l,i,j),  \<1<L' 

^  Z'  +  1</<Z 


te/),  !./<!' 

T'  +  1</<Z 

The  weight  function  introduced  in  (2.21)  is  taken  as 


1-p 

Pl=-~,  \<1<L' 


Pi=-^ - ,  Z'  +  1</<Z'  +  Z"  =  Z 

ri^Z" 


and  the  meaning  of  5  remains  the  same  (the  rms  error  of  WSBI  leading  edge). 

6.  A  corrected  ionospheric  model  is  formed  by  replacing 

u  — ^  u  +  6m 
n  =  Nu 
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7.  Ray  tracing  is  used  to  calculate  the  new  theoretical  leading  edge  g{f,  (p)  and  the 
optimum  elevation  angle  (cp,/ )  functions  for  all  values  of  frequency  and  azimuth  at 
which  the  experimental  leading  edge  points  are  specified. 

8.  The  generalized  discrepancy  is  calculated; 

P  =  (1  -  P)||l'  -  dl  +  ^^||«(0,0,^)  -  WoWlI  - 
Tl 

Iterations  are  stopped  if  p  and  5u  are  small  enough  (meaning  that  the  iterative  solution 
has  approached  the  true  one) .  Otherwise  the  process  returns  to  the  step  4  above. 
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2,  Numerical  Solution  of  Two  Dimensional  Fredholm  Integral  Equation  of  the  First 

Kind 


2.1.  General  Relationships 

Numerical  methods  for  solution  of  one  dimensional  (the  unknown  function  is  a 
function  of  one  variable)  integral  equation  of  the  first  kind  were  described  in  literature 
before  [Tikhonov  and  Arsenin,  1977].  My  realization  of  the  regularization  method  for 
two  dimensional  (the  unknown  function  is  a  fimction  of  two  variables)  integral  equation 
of  the  first  kind  is  described  here. 

Consider  the  integral  equation 

f(x,p)=  ^^K(x,p,s,t)y{s,t)dsdt,  (1) 

where  fix,p)  and  the  kernel  K(x,p,s,f)  are  considered  to  be  known  functions.  Introduce 
the  regularizing  functional 

+g,[f]  wdsdt,  (2) 

where  qs{s,t),  qf{s,t),and  w(5,t)  are  specified  positive  weight  functions.  Then  the 
smoothing  functional  is 
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Here  a  is  a  positive  constant  called  the  regularization  parameter.  In  order  to  find  the  y 

that  minimizes  this  functional  I  need  to  calculate  the  first  variation  of  the  functional.  Let 

us  first  find  variation  of  the  regularizer: 

rrf  dy  d  dy  d  ~\ 

6fi  =  2jJU5y  +  «,--5y  +  ?,--6y  ^dsdt 


In  evaluating  this  expression  the  divergence  theorem  was  used. 
Now 


JWOJ  K{x,  p,s',t  ')y(s',t')ds'dt'  -  / (x,  jo) p, s,  t)dxdp 


The  above  expression  must  be  zero  at  arbitrary  by .  Consequently  the  solution  y^  (s,t) 


that  minimizes  (3)  at  a  given  a  may  be  found  from  the  following  integro-differential 
equation  (Euler  equation): 


jjG{s,t,s',t')y(s',t')ds'dt' 

["5  d  d  d 

+  a  —  —  y{s,t)  =  F{s,t) 


The  equation  must  be  solved  with  the  boundary  condition 
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(5) 


dy  dy  ds 

=  0  (along  the  boundary) 

C/ij  \JL  C41 


G(s,t,s',t')=  ^^K{x,p,s,t)K{x,p,s'  ,t')dxdp  (6) 

F{s,t)=  ^^K{x,p,s,t)f{x,p)dxdp  (7) 

Let  5  be  the  a  priori  root-mean-square  uncertainty  in  measuring  f(x,p^,  then 
the  following  is  the  regularized  formulation  of  the  inverse  problem  (1): 

dxdp<?>^  ^^dxdp  (8) 


Q[y] 


->  mm 


Solution  of  this  problem  is  designed  as  follows: 

If  the  inequality 

dxdp  ^  (10) 

holds,  then  >>(5, 0=0  is  the  solution  of  (8),  (9).  Otherwise  the  integro-differential  equation 
(4)  is  solved  with  boundary  conditions  (5)  for  a  set  of  a  so  that  a  set  of  solutions 
y^  (s,t)  is  obtained.  These  solutions  are  then  used  to  find  the  so  called  generalized 
discrepancy  function 

p(a)  =  dxdp-b^  ^^dxdp  (11) 

The  optimum  regularization  parameter  is  the  root  of  the  equation 

p(a^/)  =  0  (12) 
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and 


y(sj)  =  yajisJ)  (13) 

is  the  solution  of  the  regularized  problem. 


2.2.  Regularization  Method  in  a  Rectangular  Region 

In  the  praetically  important  case  of  rectangular  region  of  integration  in  (1),  that  is  when 

the  region  of  integration  is  determined  as 

a  <  s  <b 
c<t  <d 

it  is  possible  to  rewrite  some  of  the  above  relationships  in  more  concrete  form.  I  will 
write  down  these  relationships  for  future  use  in  the  description  of  the  numerical  method. 
The  range  of  variables  x,p  remains  arbitrary. 

The  source  integral  equation  (1)  now  is 

b  d 

f(x,p)  =  jjK(x,p.s,  t^(s,tyisdt  (14) 

a  c 

The  Euler  equation  (4)  is  reformulated  as 


h  d 


I  ^G{s,t,s',t')y(s',t'yis'dt' 


d  d  d  d 

+  a  1  -  —  t)q,  (s,  Hs,  t)q,  {s,t)  — 


(15) 


y{s,t)  =  F{s,t) 


and  its  boundary  conditions  (in  accordance  with  (5))  are: 
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(16) 


ay 

1 

1  ^ 

s=a 

dy 

=  0,  c<t<d 
=  0,  a  <  s<b 


Finally  the  discrepancy  function  (1 1)  is  rewritten  as 


b  d 


n2 


P(a)  =  I  \\k{x,p,s,  t)yis,t)dsdt  dxdp-6^ 


(17) 


2.3.  Numerical  Realization  of  the  Regularization  Method  for  Two-Dimensional 

Integral  Equation  of  the  First  Kind 

In  order  to  accomplish  the  numerical  solution,  a  finite  element  approximation 
must  be  introduced. 


Si+2 

Si+1 


Si 


Si-1 


Figure  1 .  Grid  of  variables  s,t 


Let  us  use  the  rectangular  grid  si,tj  (Fig.  1)  subject  to  the  following  conventions; 
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.y,  =a,  Sf^  =  b 

t\  —  c,  t  =  d 

\<i<M 
l<j<N 

It  is  possible  to  conveniently  incorporate  the  boundary  conditions  by  introducing  four 
fictitious  grid  lines  numbered  by  i=-\,  i=M+\,  j=-\,  and j=N+\.  The  lines  are  placed 
infinitesimally  close  to  the  boundary  but  lie  outside  the  boundary  (Figure  2  illustrates 
how  the  grid  lines  /=-l  and  j— 1  are  introduced). 


tn=C-0  .  ,  t,=C 


s,=a 


'  s5=a-0 


Figure  2.  Auxiliary  grid  lines  (the  dashed  lines)  for 
incorporation  of  boundary  conditions 


Then  the  boundary  conditions  (16)  are  equivalent  to  the  relationships: 

y{M  +  \J)=yiM,j),  \<j<N 
t(^,0)  =  y{i,l) 

yii, N  +  l)  =  y(i, N),  \<i<  M 

All  integrals  over  are  approximated  by  finite  sums  as  follows: 
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(18) 


b  d 

a  c  \<i<M 

i^j<N 

where  r^  is  a  weight  function  which  is  specific  for  each  method  of  numerical  integration. 
1  am  using  the  trapezoid  formula,  so  that 

'•!/=^(''■+''MX^■+Vl)  (19) 

where 

=  ^'z+i  -  Si,  l<i<M 

Note  that  effectively  in  application  to  the  integration  formulas  (18)  and  (19)  there 
^0  =  =  ^0  -  -  0  . 

Now  consider  approximation  of  integrals  over  the  variables  x,p.  It  is  more 
convenient  to  refrain  from  specifying  any  structural  detail  of  the  grid  for  these  variables. 
Suppose  that  nodes  of  this  grid  are  numbered  by  one  integer  /, 

\<1<L 

then 

^Pi)Pl  ■  (21) 

/=! 

Here  Pi  is  an  appropriate  weight  function.  In  order  to  simplify  notations  in  the 
discussion  of  numerical  method  that  follows  below  I  will  replace  the  combination  Xi,pi 
by/: 

Xi,p/  ->  /. 
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I  will  use  the  following  replacements  as  well 


s,-  i 

^  j 

Here  are  the  finite  dimensional  approximations  of  (5),  (6),  and  (15): 

L 

G(i,J.i'J')  =  '^K(l,ljyC^,i',j')Pi  (22) 

/=] 

L 

=  (23) 

/=1 


J')yO ' ,  J')n'f  +  aw(i,j)y(i,  j) 


I  j 


2a  lq^(i  +  \/2,j)w(i  +  l/2,j)^ 


+  K-\ 


h 


[y(/  +  l,7)-Kh7)] 


q^(i-\l2,j)w{i-\l2,j)r  . 


h. 


7-1 


2a  g/  {Uj  + 1  /  2)w(i,y  + 1  /  2)  p  . 


-,,(,-,y-i/2)»-(,-.y-i/2)^^^.  ._i)l 


^y-1 


(24) 


It  is  advisable  to  multiply  (24)  by  /  Tq  ,  where  Tq  is  a  normalizing  constant.  This 

operation  will  transform  (24)  into  a  system  of  equations  with  symmetric  matrix  (the 
symmetry  is  with  respect  to  swapping  of  i,j  with  i',j' ,  note  that  matrix  G  is  symmetric 
in  this  sense).  The  preservation  of  the  symmetry  of  the  matrix  is  essential  for  the  sake  of 
economizing  the  use  of  computer  memory.  Finally  (24)  may  be  rewritten  as 
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(25) 


M  N 

•yQ'J')  =  F{i,j)rij  / ro  ,  \<i<M-,  1  < 7  < N 

/'=i  j'=\ 

Where 


Matrix  C  may  be  represented  in  terms  of  three 
CiiJr  =  Cff5,-,,5  Jf.  +  C?8,,.,|8,,.  +  c|6,,,6j,,„ 


M  by  N  matrices: 

+  C,y  5,-^.]  ,-,8 Jj,  +  C,y  5,y6 y+jy 


(26) 


Where 


1  ^+^-1 
2  ro/2,._, 


2  <i<M 
l<j<N 


(27) 


1  +  hi-\ 

2  ^oAy-i 


1 

^[hj-T\qt 


at 


1  <  i  <  M 
2<j<N 


(28) 


c!i=K,jwOJ)-Cl-Cl-Cfj,, 

2  3 

Cjj  and  Cij  are  zero  outside  the  range  of  indexes  indicated  in  (27),  (28).  A  possible 
choice  of  rg  is: 

(b  -a)  (d  - c) 

'  °  ir~ 

Thus,  (25)  establishes  a  system  oi  Mx  N  linear  equations  on  Mx  N  unknowns 
y{i,j) .  This  equation  is  solved  numerically  using  standard  procedure  for  systems  of 
equations  with  symmetric  matrix.  The  solution  obtained  is  an  approximate  solution  of  the 
equation  (15).  A  regularized  solution  is  obtained  as  it  is  described  in  Sections  1  and  2 
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with  evident  replacements  of  the  integrals  by  their  finite  sum  approximations  (18)-(21) 
and  equation  (15)  by  equation  (25). 
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3.  Reconstruction  of  the  Ionosphere  from  WSBI  Data  Collected  hy  the  Virginia 

OTHR 

The  Navy  OTHR  systems  are  comprised  of  the  surveillance  radar  itself,  a  quasi¬ 
vertical  sounder,  and  a  backscatter  sounder.  It  provides  nearly  simultaneous 
measurements  of  the  WSBI  and  QVI  ionograms.  The  sounder  obtains  WSBIs  for  eight 
azimuthal  beam  directions  each  separated  by  8°  (this  was  changed  to  10°  in  November 
1995).  The  central  beam  direction  of  the  sounder,  or  boresight,  is  at  175°,  i.e.  almost  due 
south. 

We  present  here  eight  examples  wherein  parameterization  (3)  was  used  in  each 
case.  This  parameterization  produces  a  model  that  has  the  same  vertical  profile  shape 
everywhere,  which  is  not  valid  for  a  real  ionosphere.  Nevertheless,  the  parameterization 
has  proven  to  be  very  useful.  It  was  demonstrated  on  simulated  and  real  experiments 
[Fridman  and  Fridman,  1994,  Fridman  et  al,  1994]  that  the  method  produces  an 
accurate  reconstruction  of  horizontal  ionospheric  variations  near  the  typical  altitude  of 
reflection  of  rays  that  form  the  leading  edge.  In  the  examples  considered  here,  the 
reflection  points  were  rather  close  to  the  F-region  peak  (lower  by  5  -  25  km)  and  for  that 
reason  we  present  our  inversions  in  terms  of  electron  density  (or  plasma  frequency)  at  the 
altitude  of  the  major  maximum  of  the  overhead  profile. 

Figures  la,  lb  contain  eight  panels  showing,  as  a  function  of  frequency  and  delay, 
the  measured  (crosses)  and  iterated  (asterisks  and  squares)  leading  edges  for  each  of  the 
azimuthal  sectors  for  backscatter  soundings  acquired  on  December  8,  1994  starting  at 
2244  UT  (UT  is  5  hours  ahead  of  local  time).  All  of  the  measured  data  points  in  these 
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panels  were  used  as  input  to  the  WSBI  inversion  technique;  the  measurement  error  for  the 
one  way  group  range  was  specified  as  6=50  km.  The  iterative  solution  starts  from  a 
model  that  assumes  a  horizontally  uniform  ionosphere  with  the  vertical  profile  being 
/7o(/z)  everywhere.  The  theoretical  leading  edge  synthesized  on  the  first  iteration  is 
shown  in  Fig.  1  by  asterisks;  the  small  squares  denote  the  leading  edge  for  the 
ionospheric  model  obtained  from  the  final  iteration. 

Figure  2  is  a  contour  plot  of  the  reconstructed  distribution  of  plasma  fi-equency  at 
the  altitude  /i=235  km  (i.e.  the  altitude  of  the  F-layer  maximum  as  obtained  by  the  QVI  at 
the  OTHR  site)  as  a  function  of  geographic  coordinates.  As  denoted  by  the  color  bar, 
isolines  on  this  and  succeeding  figures  correspond  to  integer  values  of  the  plasma 
frequency  (in  MHz);  the  red  contours  show  continental  and  island  coastlines.  An 
azimuthal  equidistant  projection  (preserving  azimuth  and  range)  from  the  OTHR  was 
used  in  creating  this  contour  map  and  quasi  Cartesian  coordinates  (in  km)  are  also 
indicated  for  convenience.  The  range  extent  of  the  area  for  which  the  reconstruction  is 
valid  is  restricted  because,  at  sufficiently  large  distances,  there  are  no  leading  edge  related 
rays  that  pass  through  the  F-region.  This  restriction  was  taken  into  account  in 
determining  the  greatest  range  extent  of  the  ionospheric  maps  presented  here. 

In  this  example,  acquired  near  18h  LT,  the  day-night  transition  pattern  is  clearly 
evident  with  a  relative  difference  of  a  factor  of  two  in  plasma  density  across  the 
terminator  region.  A  similar  example  for  data  acquired  in  March  1996  (Fig.  11)  exhibits 
approximately  the  same  difference  in  density  between  the  sunlit  and  nightside 
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3-D  Inversion  of  WSBI  Leading  Edge 
Virginia,  December  8, 1994,  22:44 
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Figure  2.  Plasma  frequency  (MHz)  at  h  =  235  km 
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ionospheres.  These  data  also  show  an  example  of  the  expanded  coverage  (80°  sector)  that 
was  begun  in  late  1 995  at  both  OTHR  facilities. 

The  leading  edge  data  shown  on  Figure  3  demonstrate  unusual  bending  down 
shape.  Inversion  of  these  WSBIs  shown  on  Figure  4  indicates  strong  latitudinal  gradient 
of  plasma  density  that  is  perhaps  a  manifestation  of  the  equatorial  anomaly  [Kelley, 
1989].  This  example  represents  data  acquired  during  the  pre-sunset  hours  in  late  summer. 
In  the  next  diagram  (Fig.  5),  which  shows  data  acquired  Ih  48m  later,  it  is  evident  that 
the  strong  gradient  had  disappeared,  with  relatively  small  variations  of  plasma  density 
that  are  typical  of  nighttime  conditions. 

Data  shown  in  Fig.  6,  7  are  characteristic  for  conditions  with  negative  gradient  of 
electron  density  typical  for  pre-midnight  time  conditions. 

The  data  shown  in  Fig.  8  were  acquired  in  the  post  noon  time  sector;  this  case  is 
of  interest  because  the  QVI  data  show  evidence  of  traveling  ionospheric  disturbances 
(TIDs)  passing  overhead,  as  each  successive  vertical  ionogram  was  drastically  different 
from  the  previous  one  (5t  =  12m).  In  general,  the  reconstructed  ionosphere  is  very 
smooth  and  the  plasma  density  varies  less  than  20%.  While  it  is  possible  that  some 
features  of  the  weak  irregular  structure  present  in  Fig.  8  are  associated  with  the  TIDs, 
further  observations  need  to  be  carried  out  in  order  to  make  a  definitive  judgment  about 
the  effect  of  TIDs  on  the  large-scale  horizontal  ionospheric  structure. 

Figure  9  demonstrates  very  weak  and  almost  homogeneous  ionosphere  observed 
several  hours  before  local  sunrise. 
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Figure  3a 
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Figure  3b 
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3-D  Inversion  of  WSBI  Leading  Edge 
Virginia,  August  5, 1994,  22:57 
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Figure  4.  Plasma  frequency  (MHz)  at  h  =  247  km 
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3-D  Inversion  of  WSBI  Leading  Edge 
Virginia,  August  6, 1994,  00:45 
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Figure  6a 
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3-D  Inversion  of  WSBI  Leading  Edge 
Virginia,  August  6, 1994,  02:21 
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Figure  7.  Plasma  frequency  (MHz)  at  h  =  243  km 


3-D  Inversion  of  WSBI  Leading  Edge 
Virginia,  December  7, 1994,  18:50 


Figiire  8.  Plasma  frequency  (MHz)  at  h  =  226  km 
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3-D  Inversion  of  WSBI  Leading  Edge 
Virginia,  December  8, 1994,  09:05 


Figure  9.  Plasma  frequency  (MHz)  at  h  =  293  km 


Figure  10  shows  ionosphere  in  local  noon  time.  The  ionosphere  is  weakly 


inhomogeneous  at  that  time. 


4.  Discussion  and  Conclusions 

We  have  demonstrated  here  the  capabilities  of  the  WSBI  inversion  method  as  a 
viable  diagnostic  tool  for  determining  the  spatial  variation  of  medium  and  large-scale 
structures  in  the  ionosphere.  It  is  particularly  valuable  that  the  method  is  able  to  produce 
two-dimensional  snapshots  of  horizontal  structure  of  the  lower  F2  region  (when  the 
parametrization  (3)  is  used).  For  the  Virginia  OTHR  these  snapshots  cover  a  64°  sector 
of  up  to  A2000  km  radius.  None  of  other  extant  experimental  techniques  (e.g.  incoherent 
scatter  and  ionospheric  tomography)  employed  for  ionospheric  measurements  are  capable 
of  providing  such  a  comprehensive  picture  of  the  horizontal  structure  of  the  ionosphere. 

The  accuracy  of  this  method  has  been  tested  using  both  simulations  and 
experimental  tests  [Fridman  and  Fridman,  1994;  Fridman  et  al.,  1994],  which  show  that 
it  provides  accurate  profiles  of  the  F  region  below  the  F2  peak.  For  other  ionospheric 
layers  the  diagnostics  may  not  be  so  effective  and  to  further  improve  this  method,  it  may 
be  necessary  to  distinguish  between  leading  edges  associated  with  the  E-  and  F-layers. 

Currently  there  are  two  OTHRs  operating  in  the  US,  one  in  Virginia  and  another 
in  Texas,  making  WSBI  and  QVI  ionospheric  soundings  every  12  minutes.  Using 
available  computational  resources,  the  WSBI  inversion  technique  described  here  is  able 
to  perform  inversions  within  the  operational  cycle  of  the  sounders.  Thus,  implementation 
of  the  WSBI  inversion  technique  at  the  existing  OTHR  sites  can  provide  real-time 
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3-D  Inversion  of  WSBI  Leading  Edge 
Virginia,  March  21, 1996,  17:45 


Figure  10.  Plasma  frequency  (MHz)  at  h  =  227km 
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3-D  Inversion  of  WSBI  Leading  Edge 
Virginia,  March  21, 1996,  23:09 
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Figure  1 1 .  Plasma  frequency  (MHz)  at  h  =  244  km 
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monitoring  of  ionospheric  irregularity  structures  over  a  very  large  geographic  area.  In 
addition  to  direct  utilization  in  the  OTHR  coordinate  registration  process,  this  unique 
information  is  of  potential  use  to  other  users  of  ionospheric  propagation  information,  as 
well  as  providing  a  valuable  input  to  the  National  Space  Weather  Program  [Wright,  1995, 
1997].  Furthermore,  routine  analysis  of  this  information  holds  the  promise  of  extending 
our  understanding  of  the  structure  and  dynamics  of  medium  and  large-scale  ionospheric 
irregularities. 

Unfortunately  this  project  was  terminated  two  months  earlier  than  planned  and  we 
were  unable  to  finish  the  work  on  WSBI  inversion  from  simultaneous  data  of  both 
backscatter  sounders.  This  development  was  stopped  on  the  debugging  stage.  If 
accomplished  it  could  provide  enhanced  BI  inversions  for  the  overlapping  coverage 
region  of  the  two  soxmders. 

Acknowledgments:  We  are  grateful  to  L.  J.  Nickisch  and  M.  Hausman  of 
Mission  Research  Corporation  for  providing  us  the  data  and  for  productive  discussions. 
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